library(knitr)
opts_chunk$set(warning = FALSE, message = FALSE, fig.path = 'figs/', dev.args = list(bg = 'transparent'), eval = T)
library(raster)
library(sf)
library(tidyverse)
library(mapview)
library(proj4shortcut)
library(plotly)
prj_geo <- geo_wgs84
prj_pro <- utm_wgs84('11s')
Map of in situ samples with chlorophyll concentration.
# get coverage of in situ points by date
insit <- read_csv('raw/CHL_2017_08_21.csv', col_names = F) %>%
rename(
site = X1,
lat = X2,
lon = X3,
chl = X4
) %>%
mutate(lat = ifelse(site == 34, 33.66954, lat))
coordinates(insit) <- ~ lon + lat
proj4string(insit) <- prj_geo
mapview(insit, zcol = 'chl', legend = T)
Overlap of in situ samples (red numbers) with the spatial extent of each image. Note that the spatial extent is rectangular and not all pixels in the extent have values.
# transform insit to projected
insit <- spTransform(insit, CRS(prj_pro))
# raster tiffs
imgs_pth <- list.files('img/GeoTiff_08_21/', pattern = '*\\.tif$', full.names = T)
# get raster extents
out <- list()
cmbs <- for(i in 1:length(imgs_pth)){
# raster aggregate, dissolve to polygon
rst_chk <- stack(imgs_pth[[i]]) %>%
.[[1]] %>%
aggregate(fact = 10, fun = mean)
values(rst_chk) <- ifelse(is.na(values(rst_chk)), NA, 1)
ply <- rst_chk %>%
rasterToPolygons(dissolve = T)
out[[i]] <- ply
}
# fortify polygons for plot
allply <- out %>%
map(., function(x){
dat <- fortify(x) %>%
mutate(fl = gsub('\\.[0-9]$', '', names(x@data)))
return(dat)
}) %>%
enframe %>%
unnest
# plot
p <- ggplot() +
geom_polygon(data = allply, aes(x = long, y = lat, group = fl), alpha = 0.5) +
geom_text(data = data.frame(insit), aes(x = lon, y = lat, label = site)) +
coord_equal()
ggplotly(p)
Extracted pixel values for each of four bands in each image where sample sites occurred within an image.
# extract values from raster
out <- list()
for(i in 1:length(imgs_pth)){
# cat(i, 'of', length(imgs_pth), '\n')
# extract raster cells by points in insit
rst <- stack(imgs_pth[[i]])
crs(rst) <- prj_pro
rst_chk <- rst %>%
raster::extract(insit, buffer = 0.5) %>%
enframe %>%
bind_cols(insit@data, .) %>%
select(-name) %>%
filter(map(.$value, is.matrix) %>% unlist) %>%
mutate(value = map(value, function(x){
sumpt <- x %>%
as.data.frame %>%
gather('img_bnd', 'din', everything()) %>%
group_by(img_bnd) %>%
summarise(
ncell = length(img_bnd),
din = mean(din)
)
return(sumpt)
})) %>%
unnest
# append to output
out[[i]] <- rst_chk
}
exts <- out %>%
do.call('rbind', .) %>%
remove_rownames() %>%
separate(img_bnd, c('img', 'bnd'), sep = '\\.') %>%
filter(!is.na(din)) # check this
exts
## # A tibble: 240 x 6
## site chl img bnd ncell din
## <int> <dbl> <chr> <chr> <int> <dbl>
## 1 31 189 IMG_170821_213834_0008 1 13 0.314
## 2 31 189 IMG_170821_213834_0008 2 13 0.166
## 3 31 189 IMG_170821_213834_0008 3 13 0.142
## 4 31 189 IMG_170821_213834_0008 4 13 0.141
## 5 31 189 IMG_170821_213840_0009 1 16 0.328
## 6 31 189 IMG_170821_213840_0009 2 16 0.175
## 7 31 189 IMG_170821_213840_0009 3 16 0.145
## 8 31 189 IMG_170821_213840_0009 4 16 0.144
## 9 31 189 IMG_170821_213846_0010 1 17 0.328
## 10 31 189 IMG_170821_213846_0010 2 17 0.170
## # ... with 230 more rows
Scatterplots of band values (din) with measured chlorophyll.
toplo <- exts %>%
group_by(site, chl, bnd) %>%
summarise(din = mean(din))
ggplot(toplo, aes(x = din, y = chl)) +
geom_point() +
facet_wrap(~bnd, ncol = 2) +
stat_smooth(method = 'lm') +
theme_bw()
Working with model selection:
library(caret)
library(MuMIn)
tomod <- exts %>%
mutate(bnd = paste0('b', bnd)) %>%
spread(bnd, din)
glb <- lm(chl ~ (b1 + b2 + b2 + b3 + b4)^3 + I(b1^2) + I(b2^2) + I(b3^2) + I(b4^2), tomod,
na.action = 'na.pass')
tmp <- dredge(glb, evaluate = F)
# createFolds(1:nrow(tomod), k = 5)
Map of in situ samples with chlorophyll concentration.
# get coverage of in situ points by date
insit <- read_csv('raw/CHL_2017_09_06.csv', col_names = F) %>%
rename(
site = X1,
lat = X2,
lon = X3,
chl = X4
) %>%
mutate(lon = ifelse(site == 39, -117.3597, lon))
coordinates(insit) <- ~ lon + lat
proj4string(insit) <- prj_geo
mapview(insit, zcol = 'chl', legend = T)
Overlap of in situ samples (red numbers) with the spatial extent of each image. Note that the spatial extent is rectangular and not all pixels in the extent have values.
# transform insit to projected
insit <- spTransform(insit, CRS(prj_pro))
# raster tiffs
imgs_pth <- list.files('img/GeoTiff_09_06/', pattern = '*\\.tif$', full.names = T)
# get raster extents
out <- list()
cmbs <- for(i in 1:length(imgs_pth)){
# raster aggregate, dissolve to polygon
rst_chk <- stack(imgs_pth[[i]]) %>%
.[[1]] %>%
aggregate(fact = 10, fun = mean)
values(rst_chk) <- ifelse(is.na(values(rst_chk)), NA, 1)
ply <- rst_chk %>%
rasterToPolygons(dissolve = T)
out[[i]] <- ply
}
# fortify polygons for plot
allply <- out %>%
map(., function(x){
dat <- fortify(x) %>%
mutate(fl = gsub('\\.[0-9]$', '', names(x@data)))
return(dat)
}) %>%
enframe %>%
unnest
# plot
p <- ggplot() +
geom_polygon(data = allply, aes(x = long, y = lat, group = fl), alpha = 0.5) +
geom_text(data = data.frame(insit), aes(x = lon, y = lat, label = site)) +
coord_equal()
ggplotly(p)
Extracted pixel values for each of four bands in each image where sample sites occurred within an image.
# extract values from raster
out <- list()
for(i in 1:length(imgs_pth)){
# cat(i, 'of', length(imgs_pth), '\n')
rst <- stack(imgs_pth[[i]])
crs(rst) <- prj_pro
rst_chk <- rst %>%
raster::extract(insit) %>%
data.frame(insit@data, .) %>%
gather('img_bnd', 'din', -site, -chl) %>%
na.omit
# out <- c(out, list(tmp))
out[[i]] <- rst_chk
}
exts <- out %>%
do.call('rbind', .) %>%
remove_rownames() %>%
separate(img_bnd, c('img', 'bnd'), sep = '\\.')
exts
## site chl img bnd din
## 1 1 62.62 IMG_170821_213858_0012 1 0.3521872
## 2 16 380.14 IMG_170821_213858_0012 1 0.3669254
## 3 17 295.34 IMG_170821_213858_0012 1 0.2611819
## 4 1 62.62 IMG_170821_213858_0012 2 0.1900249
## 5 16 380.14 IMG_170821_213858_0012 2 0.1846308
## 6 17 295.34 IMG_170821_213858_0012 2 0.1613828
## 7 1 62.62 IMG_170821_213858_0012 3 0.1485732
## 8 16 380.14 IMG_170821_213858_0012 3 0.1278632
## 9 17 295.34 IMG_170821_213858_0012 3 0.1200351
## 10 1 62.62 IMG_170821_213858_0012 4 0.1408067
## 11 16 380.14 IMG_170821_213858_0012 4 0.1263729
## 12 17 295.34 IMG_170821_213858_0012 4 0.1322991
## 13 1 62.62 IMG_170906_182425_0023 1 0.3257019
## 14 16 380.14 IMG_170906_182425_0023 1 0.3136333
## 15 1 62.62 IMG_170906_182425_0023 2 0.1682546
## 16 16 380.14 IMG_170906_182425_0023 2 0.1731694
## 17 1 62.62 IMG_170906_182425_0023 3 0.2130444
## 18 16 380.14 IMG_170906_182425_0023 3 0.2025952
## 19 1 62.62 IMG_170906_182425_0023 4 0.1939986
## 20 16 380.14 IMG_170906_182425_0023 4 0.1948591
## 21 1 62.62 IMG_170906_182431_0024 1 0.3138918
## 22 16 380.14 IMG_170906_182431_0024 1 0.3464889
## 23 1 62.62 IMG_170906_182431_0024 2 0.1581172
## 24 16 380.14 IMG_170906_182431_0024 2 0.1802249
## 25 1 62.62 IMG_170906_182431_0024 3 0.2159136
## 26 16 380.14 IMG_170906_182431_0024 3 0.2001033
## 27 1 62.62 IMG_170906_182431_0024 4 0.2123534
## 28 16 380.14 IMG_170906_182431_0024 4 0.1866015
## 29 1 62.62 IMG_170906_182437_0025 1 0.3141099
## 30 16 380.14 IMG_170906_182437_0025 1 0.2983865
## 31 31 254.50 IMG_170906_182437_0025 1 0.3381861
## 32 1 62.62 IMG_170906_182437_0025 2 0.1613249
## 33 16 380.14 IMG_170906_182437_0025 2 0.1533507
## 34 31 254.50 IMG_170906_182437_0025 2 0.1642923
## 35 1 62.62 IMG_170906_182437_0025 3 0.2072103
## 36 16 380.14 IMG_170906_182437_0025 3 0.1837897
## 37 31 254.50 IMG_170906_182437_0025 3 0.2112308
## 38 1 62.62 IMG_170906_182437_0025 4 0.1993446
## 39 16 380.14 IMG_170906_182437_0025 4 0.1752855
## 40 31 254.50 IMG_170906_182437_0025 4 0.2066646
## 41 1 62.62 IMG_170906_182443_0026 1 0.3077186
## 42 16 380.14 IMG_170906_182443_0026 1 0.3145586
## 43 31 254.50 IMG_170906_182443_0026 1 0.3824544
## 44 1 62.62 IMG_170906_182443_0026 2 0.1680164
## 45 16 380.14 IMG_170906_182443_0026 2 0.1602393
## 46 31 254.50 IMG_170906_182443_0026 2 0.1722880
## 47 1 62.62 IMG_170906_182443_0026 3 0.2157355
## 48 16 380.14 IMG_170906_182443_0026 3 0.2072434
## 49 31 254.50 IMG_170906_182443_0026 3 0.2396860
## 50 1 62.62 IMG_170906_182443_0026 4 0.2234698
## 51 16 380.14 IMG_170906_182443_0026 4 0.1862677
## 52 31 254.50 IMG_170906_182443_0026 4 0.2482327
## 53 1 62.62 IMG_170906_182449_0027 1 0.3003604
## 54 16 380.14 IMG_170906_182449_0027 1 0.3132594
## 55 17 295.34 IMG_170906_182449_0027 1 0.3272939
## 56 31 254.50 IMG_170906_182449_0027 1 0.4016975
## 57 1 62.62 IMG_170906_182449_0027 2 0.1576792
## 58 16 380.14 IMG_170906_182449_0027 2 0.1483788
## 59 17 295.34 IMG_170906_182449_0027 2 0.1645345
## 60 31 254.50 IMG_170906_182449_0027 2 0.1726661
## 61 1 62.62 IMG_170906_182449_0027 3 0.2268322
## 62 16 380.14 IMG_170906_182449_0027 3 0.2526465
## 63 17 295.34 IMG_170906_182449_0027 3 0.2137336
## 64 31 254.50 IMG_170906_182449_0027 3 0.2233939
## 65 1 62.62 IMG_170906_182449_0027 4 0.1983301
## 66 16 380.14 IMG_170906_182449_0027 4 0.2199366
## 67 17 295.34 IMG_170906_182449_0027 4 0.2050995
## 68 31 254.50 IMG_170906_182449_0027 4 0.2229731
## 69 17 295.34 IMG_170906_182455_0028 1 0.2903806
## 70 31 254.50 IMG_170906_182455_0028 1 0.4145821
## 71 32 214.83 IMG_170906_182455_0028 1 0.0000000
## 72 17 295.34 IMG_170906_182455_0028 2 0.1540359
## 73 31 254.50 IMG_170906_182455_0028 2 0.1831722
## 74 32 214.83 IMG_170906_182455_0028 2 0.0000000
## 75 17 295.34 IMG_170906_182455_0028 3 0.1966201
## 76 31 254.50 IMG_170906_182455_0028 3 0.2158420
## 77 32 214.83 IMG_170906_182455_0028 3 0.2106364
## 78 17 295.34 IMG_170906_182455_0028 4 0.2066191
## 79 31 254.50 IMG_170906_182455_0028 4 0.2331240
## 80 32 214.83 IMG_170906_182455_0028 4 0.1957703
## 81 2 486.96 IMG_170906_182501_0029 1 0.3103488
## 82 17 295.34 IMG_170906_182501_0029 1 0.2969343
## 83 32 214.83 IMG_170906_182501_0029 1 0.4165199
## 84 2 486.96 IMG_170906_182501_0029 2 0.1686566
## 85 17 295.34 IMG_170906_182501_0029 2 0.1491194
## 86 32 214.83 IMG_170906_182501_0029 2 0.1849540
## 87 2 486.96 IMG_170906_182501_0029 3 0.2082499
## 88 17 295.34 IMG_170906_182501_0029 3 0.2056280
## 89 32 214.83 IMG_170906_182501_0029 3 0.2059632
## 90 2 486.96 IMG_170906_182501_0029 4 0.2313217
## 91 17 295.34 IMG_170906_182501_0029 4 0.2107968
## 92 32 214.83 IMG_170906_182501_0029 4 0.2166646
## 93 2 486.96 IMG_170906_182507_0030 1 0.3085389
## 94 17 295.34 IMG_170906_182507_0030 1 0.3139217
## 95 18 305.90 IMG_170906_182507_0030 1 0.2996119
## 96 32 214.83 IMG_170906_182507_0030 1 0.4486934
## 97 2 486.96 IMG_170906_182507_0030 2 0.1510785
## 98 17 295.34 IMG_170906_182507_0030 2 0.1526142
## 99 18 305.90 IMG_170906_182507_0030 2 0.1638244
## 100 32 214.83 IMG_170906_182507_0030 2 0.1995005
## 101 2 486.96 IMG_170906_182507_0030 3 0.1965176
## 102 17 295.34 IMG_170906_182507_0030 3 0.2199249
## 103 18 305.90 IMG_170906_182507_0030 3 0.2149738
## 104 32 214.83 IMG_170906_182507_0030 3 0.2212107
## 105 2 486.96 IMG_170906_182507_0030 4 0.2272809
## 106 17 295.34 IMG_170906_182507_0030 4 0.2050937
## 107 18 305.90 IMG_170906_182507_0030 4 0.1843024
## 108 32 214.83 IMG_170906_182507_0030 4 0.2479621
## 109 2 486.96 IMG_170906_182513_0031 1 0.3167584
## 110 18 305.90 IMG_170906_182513_0031 1 0.3033060
## 111 32 214.83 IMG_170906_182513_0031 1 0.4871817
## 112 2 486.96 IMG_170906_182513_0031 2 0.1560977
## 113 18 305.90 IMG_170906_182513_0031 2 0.1605711
## 114 32 214.83 IMG_170906_182513_0031 2 0.1999720
## 115 2 486.96 IMG_170906_182513_0031 3 0.2126312
## 116 18 305.90 IMG_170906_182513_0031 3 0.2118794
## 117 32 214.83 IMG_170906_182513_0031 3 0.2445490
## 118 2 486.96 IMG_170906_182513_0031 4 0.2171975
## 119 18 305.90 IMG_170906_182513_0031 4 0.2094065
## 120 32 214.83 IMG_170906_182513_0031 4 0.2229194
## 121 2 486.96 IMG_170906_182519_0032 1 0.3505707
## 122 18 305.90 IMG_170906_182519_0032 1 0.3049344
## 123 33 151.64 IMG_170906_182519_0032 1 0.4672850
## 124 2 486.96 IMG_170906_182519_0032 2 0.1684053
## 125 18 305.90 IMG_170906_182519_0032 2 0.1514488
## 126 33 151.64 IMG_170906_182519_0032 2 0.2076212
## 127 2 486.96 IMG_170906_182519_0032 3 0.2052642
## 128 18 305.90 IMG_170906_182519_0032 3 0.2016440
## 129 33 151.64 IMG_170906_182519_0032 3 0.2221849
## 130 2 486.96 IMG_170906_182519_0032 4 0.2000579
## 131 18 305.90 IMG_170906_182519_0032 4 0.1999472
## 132 33 151.64 IMG_170906_182519_0032 4 0.2175644
## 133 18 305.90 IMG_170906_182525_0033 1 0.3112594
## 134 33 151.64 IMG_170906_182525_0033 1 0.5144554
## 135 18 305.90 IMG_170906_182525_0033 2 0.1594713
## 136 33 151.64 IMG_170906_182525_0033 2 0.2494634
## 137 18 305.90 IMG_170906_182525_0033 3 0.2018453
## 138 33 151.64 IMG_170906_182525_0033 3 0.2550449
## 139 18 305.90 IMG_170906_182525_0033 4 0.2235672
## 140 33 151.64 IMG_170906_182525_0033 4 0.2828100
## 141 3 184.99 IMG_170906_182531_0034 1 0.3149959
## 142 19 292.94 IMG_170906_182531_0034 1 0.3125520
## 143 33 151.64 IMG_170906_182531_0034 1 0.5087944
## 144 3 184.99 IMG_170906_182531_0034 2 0.1513753
## 145 19 292.94 IMG_170906_182531_0034 2 0.1535808
## 146 33 151.64 IMG_170906_182531_0034 2 0.2157703
## 147 3 184.99 IMG_170906_182531_0034 3 0.1921439
## 148 19 292.94 IMG_170906_182531_0034 3 0.2061476
## 149 33 151.64 IMG_170906_182531_0034 3 0.2267226
## 150 3 184.99 IMG_170906_182531_0034 4 0.1936413
## 151 19 292.94 IMG_170906_182531_0034 4 0.2079446
## 152 33 151.64 IMG_170906_182531_0034 4 0.2290599
Scatterplots of band values (din) with measured chlorophyll.
toplo <- exts %>%
group_by(site, chl, bnd) %>%
summarise(din = mean(din))
ggplot(toplo, aes(x = din, y = chl)) +
geom_point() +
facet_wrap(~bnd, ncol = 2) +
stat_smooth(method = 'lm') +
theme_bw()